home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / TQLI.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  76 lines
  1. PROCEDURE tqli(VAR d,e: glnp; n: integer; VAR z: glnpnp);
  2. (* Programs using routine TQLI must define the types
  3. TYPE
  4.    glnp = ARRAY [1..np] OF real;
  5.    glnpnp = ARRAY [1..np,1..np] OF real;
  6. where np is the physical dimension of the matrix to be analyzed. *)
  7. LABEL 1,2;
  8. VAR
  9.    m,l,iter,i,k: integer;
  10.    s,r,p,g,f,dd,c,b: real;
  11. FUNCTION sign(a,b: real): real;
  12.    BEGIN
  13.       IF (b < 0) THEN sign := -abs(a) ELSE sign := abs(a)
  14.    END;
  15. BEGIN
  16.    IF  (n > 1)  THEN BEGIN
  17.       FOR i := 2 TO n DO BEGIN
  18.          e[i-1] := e[i]
  19.       END;
  20.       e[n] := 0.0;
  21.       FOR l := 1 TO n DO BEGIN
  22.          iter := 0;
  23. 1:         FOR m := l TO n-1 DO BEGIN
  24.             dd := abs(d[m])+abs(d[m+1]);
  25.             IF  (abs(e[m])+dd = dd) THEN  GOTO 2
  26.          END;
  27.          m := n;
  28. 2:         IF (m <> l) THEN BEGIN
  29.             IF (iter = 30) THEN BEGIN
  30.                writeln('pause in routine TQLI');
  31.                writeln('too many iterations'); readln
  32.             END;
  33.             iter := iter+1;
  34.             g := (d[l+1]-d[l])/(2.0*e[l]);
  35.             r := sqrt(sqr(g)+1.0);
  36.             g := d[m]-d[l]+e[l]/(g+sign(r,g));
  37.             s := 1.0;
  38.             c := 1.0;
  39.             p := 0.0;
  40.             FOR i := m-1 DOWNTO l DO BEGIN
  41.                f := s*e[i];
  42.                b := c*e[i];
  43.                IF (abs(f) >= abs(g)) THEN BEGIN
  44.                   c := g/f;
  45.                   r := sqrt(sqr(c)+1.0);
  46.                   e[i+1] := f*r;
  47.                   s := 1.0/r;
  48.                   c := c*s
  49.                END ELSE BEGIN
  50.                   s := f/g;
  51.                   r := sqrt(sqr(s)+1.0);
  52.                   e[i+1] := g*r;
  53.                   c := 1.0/r;
  54.                   s := s*c
  55.                END;
  56.                g := d[i+1]-p;
  57.                r := (d[i]-g)*s+2.0*c*b;
  58.                p := s*r;
  59.                d[i+1] := g+p;
  60.                g := c*r-b;
  61.             (* Next loop can be omitted if eigenvectors not wanted *)
  62.                FOR k := 1 TO n DO BEGIN
  63.                   f := z[k,i+1];
  64.                   z[k,i+1] := s*z[k,i]+c*f;
  65.                   z[k,i] := c*z[k,i]-s*f
  66.                END
  67.             END;
  68.             d[l] := d[l]-p;
  69.             e[l] := g;
  70.             e[m] := 0.0;
  71.             GOTO 1
  72.          END
  73.       END
  74.    END
  75. END;
  76.